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1.  INTRODUCTION 


In  two  previous  reports,  Simkins  et  al  [1],  [2]  derived  the 
equation  of  motion  of  a simply  supported,  thin  elastic  beam  subjected 
to  an  internal  pressure  (sometimes  called  a Bourdon  force)  and  a 
concentrated  mass  moving  at  uniform  velocity  along  its  axis.  The 
forced  dynamic  response  of  the  system  was  determined  with  the  aid  of 
the  finite  element  method.  It  was  pointed  cut  in  these  reports  that 
if  the  moving  mass  is  an  appreciable  fraction  of  the  total  mass  of 
the  structure  or  if  the  velocity  of  the  mass  approaches  a critical 
value,  the  dynamic  deflections  of  the  system  can  become  significantly 
large.  In  addition,  if  the  moving  mass  remains  on  the  beam  for  a 
sufficiently  long  period  of  time,  the  question  of  the  stability  of 
the  system  must  be  examined. 

It  is  the  objective  of  this  report,  then,  to  determine  the 
fundamental  zone  of  instability  in  the  mass-velocity  plane.  Because 
of  the  complexity  of  the  partial  differential  equation  of  motion,  one 
cannot  hope  to  determine  easily  the  exact  solution  of  the  problem. 
However,  through  use  of  the  Galerkin  procedure,  it  is  possible  to 
obtain  a useful  approximate  solution  of  the  problem. 


2.  THE  EQUATION  OF  MOTION 

The  partial  differential  equation  of  motion  of  a slender  beam  of 
length  l,  uniform  cross  sectional  area  A,  and  mass  density  p subjected 
to  a gravity  body  force  and  a concentrated  mass  m moving  at  uniform 
velocity  along  its  axis  (the  x-j  - axis)  is 
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EIa)  im  + t f + 2Vu  + V2(j  + g)  6(x,  - Vt, ) + 
•111!  ’It]  ’ll  I 1 

+ pAw  = pgA,  o<x1  <£.,  (1) 

’*1*1  1 

where  u>  = 3oj/3x.  w . = 3u)/3ti  (t,  denotes  time),  El  is  the 

>1  I > * U«|  1 I 

flexural  rigidity  of  the  beam,  g the  acceleration  of  gravity,  and 
6(xj  - Vt,)  the  Dirac  delta  function.  Since  the  beam  is  simply 
supported  at  both  ends,  the  apposite  boundary  conditions  are 


uo  = (jo  n _ = o at  x,  = o £, 

>11  I > 


(2) 


The  coefficient  of  the  Dirac  delta  function  in  equation  (1)  represents 
the  influence  of  the  concentrated  mass  m moving  along  the  x-j-axis. 

The  objective  of  the  analysis  to  follow  consists  of  determining  the 
values  of  the  mass  m and  velocity  V parameters  for  which  the  system 
becomes  unstable.  The  term  representing  the  effect  of  a uniform 
internal  pressure  has  not  been  included  in  equation  (1).  Some  aspects 
of  the  boundary  value  problem  stated  above  have  been  discussed  by 
Bolotin  [3]  and  Inglis  [4].  An  interesting  related  problem  has  been 
solved  by  Stanisic  and  Hardin  [5].  Fryba  [6]  has  published  a book  that 
summarized  the  principal  results  of  a fairly  broad  class  of  moving 
load  problems. 

If  the  changes  of  variables 

“(xi  ,ti ) ~ £y(xft)  x-j  = Six  t-|  = at 

are  made  and  the  parameters 

a2  = pAJlVEI  p = m/pAi  M = PgAH3/EI  v = Va/l  (3) 

» » > > 
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are  introduced,  then  the  boundary  value  problem  in  equations  (1) 
and  (2)  can  be  expressed  in  dimensionless  form  as 

y,v  + p(M  + y + 2vy'  + v2y")  6(x  - vt)  + y - M = o o<x<l  (4) 

> 9 

y = y"  = o at  x = o 1 (5) 

9 9 

where  y'  = 3y/3x,  y = 3y/3t,  etc. 

3 . THE  GALERKI N PROCEDURE 

In  view  of  the  boundary  conditions  in  equation  (5),  one  is 
motivated  to  assume  a solution  for  equation  (4)  in  the  form 

oo 

y(x,t)  = l qn(t)  sin(nTTx),  (6) 

n=l 

where  the  generalized  coordinates  qn  are  unkown  functions  of  time. 
Substitution  of  equation  (6)  into  equation  (4),  multiplication  of  the 
t result  by  sin(kirx),  and  integration  of  the  product  yield 

N 

q,  + ( krr ) **  q.  + 2p  7 [q  sin  n0  + 2mrvq  cos  n0  - 
k k K ^ nn  n 

- (mrv)2  q sin  ne]  sin  ke  = ~ [1  - (-1 ) ] - 
n ku 

i»  - 2pM  sin  k0,  (7) 


in  which  q.  = d2q,/dt2  and  0 = irvt. 


which  simply  represents  a superposition  of  the  complementary  and  the 
particular  solutions.  For  the  purpose  of  determining  the  state  of 
stability  of  the  system  under  consideration,  it  is  sufficient  to 
examine  the  complementary  solution  only  which  satisfies  the  ordinary 
differential  equation 

N 

q + ( k-rr) **  q + 2p  7 [q  sin  n6  + 2nTrvq  cos  n0  - 
k k n^i  n n 

- (nirv)2  q sin  ne]  sin  k0  = 0,  (8) 

n 

where  the  superscript  c has  been  deleted. 

If  N = 1 , equation  (8)  becomes  simply 

^1  + Tt^q-j  + p[q-j(l  - cos  20)  + 2tTvq^  sin  20  - 

- ( ttv ) 2q^  (1  - cos  20)]  = 0.  (9) 

For  N = 2,  one  has  a pair  of  coupled  differential  equations,  namely, 
q^  + 7r4q-|  + p[q^(l  - cos  20)  + 2 vq^  sin  20  - 

- (iTv)2q1  (1  - cos  20)  + q2  (cos  0 - cos  30)  + 

+ 47ivq2  (sin  30  - sin0)  - (2iTv)2q2  (cos  0 - cos  30)]  = 0 (10) 

and 

q^  + (27r)4q2  + pfq^  (cos  0-  cos  30)  + 2™^  (sin  0 + sin  30)  - 

- (ttv ) 2q ^ (cos  0 - cos  30)  + q2  (1  - cos  40 ) + 

+ 4uvq2  sin  40  - (2fTv)2q2  (1  - cos  40)]  = 0.  (11) 
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4.  THE  STABILITY  BOUNDARIES 
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A useful  approximation  for  the  stability  boundaries  in  the  vp- 
plane  can  be  obtained  from  the  solution  of  equation  (9).  If  more 
refined  results  might  be  required,  one  must  solve  equations  (10)  and 
(11),  to  which  the  methods  of  references  [7]  and  [8]  probably  may  be 
applied. 


If  one  defines  w = 7iv,  0 = cot,  y = 7iz/co,  and  q^(t)  = q(0),  then 


equation  (9)  can  be  expressed  as 

<2"  — ~ rd2_9.  (1  - cos  2 ) + sin  2 - q ( 1 - cos  2 ) 

do2  d0 


d-4  + Yzq  + u 
d02 


= 0,  (12) 


which  is  related  to  a Mathieu  or  Hill  type  of  ordinary  differential 
equation  with  periodic  coefficients.  To  solve  equation  (12),  it  is 
efficacious  to  use  Whittaker's  technique  [9]  which  was  first  used  for 
the  purpose  of  solving  Mathieu's  equation.  This  technique  exploits 
the  fact  that  equation  (12)  has  a solution  of  the  form 

q(0)  = erG  4>(0),  (13) 

where  $(0)  is  a periodic  function.  Therefore,  the  state  of  stability 
of  the  system  can  be  ascertained  once  the  coefficient  r in  the 
exponential  factor  of  the  solution  in  equation  (13)  is  determined  as 
a function  of  the  parameters  of  the  system.  In  particular,  under  the 
assumption  that  r is  real,  the  motion  consists  of  (i)  oscillations 
with  an  exponentially  decaying  amplitude  (asymptotic  stability)  if 
r<o,  (ii)  oscillations  of  constant  amplitude  (neutral  stability)  if 
r=o,  and  (iii)  oscillations  with  exponentially  growing  amplitude 
(instability)  if  r>o.  Thus,  the  boundaries  of  regions  of  stability 
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in  stability  maps  must  be  obtained  from  the  condition 


r = 0.  (14) 

Substitution  of  equation  (13)  into  equation  (12)  yields 

<{)"  + 2r<t>‘  + (r2  + y2 )<J>  + y[(4>"  + 2r4>'  + r2<J>)(l  - cos  20)  + 

+ 2(r<j>  + <}>’)  sin  20  - (1  - cos  20)]  = 0,  (15) 

■where  <J>'  = dcp/de , etc.  If  the  mass  ratio  y,  as  defined  in  equation 
(3),  is  a small  number  relative  to  unity,  a regular  perturbation 
procedure,  employing  y as  the  expansion  parameter,  can  be  applied  to 
equation  (15).  Of  course,  y<l  whenever  the  mass  of  the  moving 
particle  is  less  than  the  mass  of  the  beam.  If  y>l , then  it  may  be 
possible  to  use  a perturbation  method  upon  introducing  a new  small 
parameter,  say,  e defined  as  e = y/(l  + y),  where  o<e<l . Such  a 
scheme  has  been  used  in  reference  [10].  Throughout  the  remainder  of 
this  report,  it  will  be  assumed  that  y<l . 

Hence,  one  may  suppose  that 

oj2  = a2[l  + yi?i  + y262  + (y3)],  (16) 

where  a,  B]  , B2>  • • • are  unknowns,  so  that 

y2  = 7T**/co2  = (n2/a)2  LI  + UB1  + y2B2  + (y3)]’1 

= Y2  [1  - VB-j  + U2(62  * B2)  + (y3)],  (17) 

with  y = nz/a.  Furthermore,  it  is  assumed  that 


r = yr]  + y2r?  + (y3), 


(19) 


cf  ( 0 ) = sin  (yG  - 0Q)  + (0 ) + p Z<J>2 ( 0 ) + (u3)» 

where  0q  is  an  arbitrary  phase  angle.  This  set  of  formulas  yields 
the  exact  expression  for  q(o)  when  p = o. 

Substitution  of  equations  (17)  to  (19)  into  equations  (15)  and 
utilization  of  the  familiar  steps  of  the  classical  regular  perturbation 
procedure,  one  obtains,  writing  ip  = yO  - 0 , 

0"  + y = - 2r^y  cos  ^ + (1  + y2  + y2 By ) sin  ^ - 

- ^(1  + y)2  sin(ip  + 20)  - 1 (1  - y)2  s i n ( 0 - 26),  (20) 

<J>2  + = “ rl  sin  ^ “ Zr2y  C0S  ^ ~2ri<I>l  + ^l'h  ~ 

- y2  (B2  - B2)  Sin  i|>  - (<j>"  + 2r^y  cos  xp) ( 1 - cos  26)  - 

- 2 (cp j + r^  sin  \p)  sin  26  + <t>^  (1  - cos  20),  (21) 

etc.  Equations  (20)  and  (21)  are  non-homogeneous  ordinary  differential 
equations  that  can  be  solved  by  elementary  methods. 

The  arguments  of  the  trigonometric  terms  on  the  right  side  of 
equation  (20)  are 

{p  + 29  = (y  + 2)6  - 6q  and  ip  - 20  = (y  - 2)0  - 0^. 

If,  in  particular,  y - 2 = - y,  then  y = 1,  which  implies  that  a = tt2, 
and 

ip  + 20  = 30  - 0Q  and  ip  - 29  = -(0  + 0 ).  (22) 
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Substitution  of  equation  (22)  into  equation  (20)  yields,  in  this 
ci rcumstance. 


d 


+ ♦1  = - 2r]  cos  (6  - e ) + (2  + ^)  sin  (9  - 0Q) 

- 2 sin  (30  - 0 ).  (23) 

To  avoid  secular  terms  in  <(>  (0),  one  must  equate  the  coefficients  of 
cos(o  - 0^)  and  sin(0  - 0 ) to  zero.  Consequently, 

r1  = 0,  = - 2,  (24) 

and 

+ $1  = - 2 sin  (30  - 0^).  (25) 

The  particular  solution  of  equation  (25)  is 

<Me)  = 7 sin  (30  - 0 ).  (26) 

i 4 U 

Substitution  of  y = 1 and  equation  (26)  into  equation  (21)  yields, 
after  some  rearrangement, 

<>2  + *2  = " (2r2  ' 8 Sin  290)  C°S  (Q  ' V + (62  " T + 

3 

+ ~ cos  20q)  sin  (0  - 0 ) + 2 sin  (30  - 0^)  - 

- 2 sin  (50  - 0Q).  (27) 


To  eliminate  secular  terms  in  $ (0),  one  sets  the  coefficients  of 
cos  (0  - 0 ) and  sin  (0  - 0 ) in  equation  (27)  equal  to  zero.  Hence, 


r 


and 


d>”  + = 2 sin  (30  - 0Q)  - 2 sin  (50  - 0Q). 

By  virtue  of  equations  (24)  and  (28),  one  has  from  equation  (18) 
r = | y2  sin  20Q  + (y3).  (29) 


Therefore,  it  follows  immediately  from  equations  (14)  and  (29)  that 

the  condition  for  neutral  stability  is  sin  26  = 0,  whence  0=0 

0 0 

and  n/2.  Hence,  by  the  second  relation  in  equation  (28), 


B0  = ^{21  - 3 cos  29  ] = /9//2 

2 4 0JeQ  = o,  tt/2  \6  if 


fV°' 

0O  = tt/2  . 


(30) 


Finally,  from  to  = ttv,  a = it2,  and  equations  (16),  (24),  and  (30),  one 
determines 


v’  = [1  * PS,  * p2S2  ♦ . 0i  „/2j 

whence  the  expressions  for  the  stability  boundaries  are  found  to  be 


V tt ( 1 - 2y  + 9y2/2) 


1/2 


and 


v 7T ( 1 - 2p  + 6y2)^^. 


(31) 

(32) 


% 

*4 


V . 
t 


hi 


Equations  (31)  and  (32)  provide  the  boundary  curves  for  the 
shaded  zone  of  instability  (i.e.,  oscillations  with  an  exponentially 
growing  amplitude)  shown  in  the  stability  map  in  the  vy-plane  in 
Figure  1.  Whenever  the  values  of  the  parameters  v and  y are  such  that 
the  point  (v/iT,y)  falls  within  this  shaded  region,  one  must  conclude 
that  the  motion  of  the  system  is  unstable.  Whenever  the  point  (v/tt,u) 


a 

f 

► 

iJL 


-J 


falls  within  the  unshaded  region,  the  motion  of  the  system  is 
asymptotical ly  stable.  Clearly,  the  region  of  instability  appears  at 
v = tt  when  y = 0.  As  the  value  of  y is  increased,  the  values  of  v 
for  instability  decrease  initially  below  v = tt,  but,  as  y becomes 
larger,  the  width  of  this  region  grows  rapidly. 

Without  determining  additional  terms  for  the  radicands  in 
equations  (31)  and  (32),  one  can  easily  derive  more  accurate  approx- 
imations to  the  sums  of  the  implied  series  than  the  original  three 
terms,  which  are  of  the  form 

(v/tt)2  = 1 - ay  + by2  + . . . , (33) 

where  a = 2 and  b = 9/2  or  6 for  equation  (31)  or  (32),  respectively. 

To  accomplish  this,  let 

S1  = 1 , = 1 - ay,  = 1 - ay  + by2 

denote  the  first  three  partial  sums  of  the  infinite  series  in 
equation  (33).  Shanks  [11]  has  shown  that  the  nqn-linear  transformatior 

(„/„)>  = hh'Jt  , -I- 

1 S1+S3-2S2  a+by 

yields  a more  accurate  value  of  the  sum  of  the  series  in  equation 
(33).  In  particular,  one  may  replace  equations  (31)  and  (32)  by 


v = ir[  (4  + y )/  (4  + 9y)]  , 


v = tt[(1  + y ) / ( 1 + ay)]'7", 


respecti vely. 


In  Figure  2,  a stability  map  that  is  more  accurate  than  the  one 
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igure  2.  Ref 


given  in  Figure  1 is  plotted.  The  boundary  curves  for  the  shaded 
zone  of  instability  are  those  obtained  from  equations  (34)  and  (35). 
Actually,  a close  comparison  of  Figures  1 and  2 reveals  that  the 
portions  of  the  stability  boundaries  obtained  from  equations  (31)  - 
(32)  and  (34)  - (35)  are  in  excellent  agreement  for  0<u<0.1.  For 
y>0.1,  however,  the  curves  in  Figure  2 do  not  show  the  marked 
reverse  curvature  that  is  evident  in  Figure  1.  Moreover,  the  width 
of  the  zone  of  instability  in  Figure  2 does  not  increase  so  rapidly 
as  that  shown  in  Figure  1. 
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